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In premixed turbulent flames the presence of intense mixing zones located in 
front of and behind the flame surface leads to a requirement to study the behavior 
of iso-concentration surfaces defined for all values of the progress variable (equal to 
unity in burnt gases and to zero in fresh mixtures). To support this study, some 
theoretical and mathematical tools devoted to level surfaces axe first developed. 
Then a database of direct numerical simulations of turbulent premixed flames is 
generated and used to investigate the internal structure of the flame brush, and a 
new pdf model based on the properties of iso-surfaces is proposed. 


1. Introduction 

In order to improve existing models of premixed turbulent combustion (Borghi 
1988), it is necessary to build a bridge from zones where intense reactive activity 
occurs in laminar flamelets to the more distributed reaction zones found in front of 
and behind these propagating flame surfaces. This is especially true if one wishes 
to predict not only the global structure of the flame, but also the concentration 
level of pollutants using any probabilistic approach. Indeed, the turbulent mixing 
acting in the preheat zone can strongly influence the structure of the entire flame; 
therefore, a description of this mixing, together with descriptions of chemical sources 
or of topological properties of the flame surface, must be used to account for their 
mutual interactions. 

The objective of this work is the development of a formalism, adapted to premixed 
turbulent combustion modeling, which brings the possibility to treat simultaneously 
flamelet behavior and out-of-flamelet mixing. This is a crucial point because they 
both interact to control the properties of the propagation velocity of the flame front 
(Bray 1990, Pocheau 1992). 

To account for the possible presence of flamelets, the flame surface density ap- 
proach (Marble and Broadwell 1977) is extended to the concept of a surface density 
function and coupled with probability density function (pdf) method (Pope 1985). 
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First some theoretical considerations and mathematical tools concerning iso- 
surfaces in fluid mechanical problems are developed. In Section 2.1 a useful math- 
ematical relation between surface and volume integrals is introduced which has 
several applications of theoretical and practical importance. It is used first (Section 
2.2) to establish the transport equation for global variables such as the total surface 
area of iso-concentration surfaces, and from this the classical expression for pre- 
mixed flame stretch is recovered. The solution of the area equation for a material 
surface is proposed (Section 2.3). In particular, it is shown that the combination 
of the introduced mathematical relation with Lagrangian variables leads to an ex- 
pression for the area of a temporally evolving level surface, which only requires the 
integration over the surface at a fixed time. In Section 2.4 the same formalism is 
used in conjunction with an integral transform for the computation of the area of 
iso-surfaces without computing the location of the surfaces. This leads to a novel 
method for the efficient computation of the surface area of iso-surfaces, which is 
superior to triangulation methods that keep track of the connectedness of the sur- 
face. Finally, the concept of a surface density function (sdf) emerges (Section 2.5) 
along with its transport equation, and a unified treatment of pdf and sdf equations 
is discussed (Section 2.6). 

A direct numerical simulation (DNS) database of premixed turbulent combustion 
is developed and used in Section 3 to study the internal structure of the turbulent 
flame, and various key quantities identified in Section 2 are analyzed. In particular, 
the behavior of quantities integrated over iso- concentration surfaces is carefully 
examined within the inner structure of the premixed flame, for surfaces located 
between fresh and burnt gases. 

The pdf formulation reduces the modeling difficulties emerging when the mean 
burning rate needs to be estimated to the finding of a closed expression for the 
small scales mixing which depends on the distribution of the gradient of the reactive 
species along their iso-surfaces. The surface density function contains information 
about these gradients, especially within the flamelet zone. Based on this remark, an 
attempt to couple the pdf with the sdf is proposed in Section 4, and the resulting 
modeled expressions are tested against their exact values extracted from the DNS. 

2. Transport equations for surfaces 

2.1 Introduction 

In multi-composition fluid mechanical problems, the starting point for studying 
the quantities characterizing the physical properties of the flow field are usually 
budget equations for volume integrals over a volume of fluid T>. To derive equations 
for the corresponding properties of iso-surfaces, it is necessary to link volume inte- 
grals and integrals over these surfaces. The following theorem of geometric measure 
theory (Maz’ja 1985, ch.1.2.4, p.37, Constantin 1992) provides the bridge between 
volume integrals and integrals over iso-surfaces. 

Theorem : Let F(x) be a Borel-measurable nonnegative function on VCR 2 3 open. 
Let *(x) e C°°{V). Then 
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OO 

f dxF(x)\V$(x)\ = f dip f dA(x)F(x) , 
v 0 S+(<p) 


( 1 ) 


where S$(<p) = {x € T> : $(x) = is the level surface or iso-surface of the 
field $(x) on which $(x) takes the value ip. This condition specifies that the level 
surfaces are contained in the domain V (no intersection of the surfaces with the 
domain boundary). It is possible to remove this restriction, but this case will not 
be pursued here. 

The co-area formula follows at once if the function F(x) is specified as: 


F(x) = f(x)6Wx)-<p) , 

where <5($(x) — p) is the diract delta function, leading to the desired result 

/ d^(x)\V^(x)\6(i(x) - <p) = / dA(x)f(x) . (2) 

T> $♦(*) 

2.2 Transport equation for global variables 

The total area of level surfaces is an example of a global variable. In general, a 
global variable is defined by 

= J dA(x)f(x, t) , (3) 

s*M 

where the unspecified function /(x, t) must satisfy the conditions of the above the- 
orem. In most cases it will be a smooth function of location and time to allow 
differentiation. Application of Eq. 2 to Eq. 3 and differentiation with respect to 
time leads to the transport equation for in terms of volume integrals 

^( V? ,<) + ^r* = 5 / + 5* . (4) 

Three contributions determine the time rate of change of the global variable 
F$(y>it): a convective transport term in scalar space and two source terms. The 
flux in scalar space is defined by 

C* = J dx6(* - V )|V#|/(*,0|£ , (5) 

V 

the source due to the evolution of /(x, t) is 

S/(vM) = j dx6($-<p) |V*|^ , 

V 


( 6 ) 
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and finally the source due to the evolution of the scalar field defining the level 
surfaces is 


S*(*M) = J dx6 ( * - . (7) 

V 

The co-area formula Eq. 2 can be applied once more to obtain an equivalent 
version in terms of surface integrals 

Tr ( * 0 - / <m (^ + 4'^ |v * I)+/ £) ■ (8) 

S*(v>) 

where is the time rate of change following points on the level surface. The time 
rate of change of the global variable is seen to be caused by the change of the 
function f{x,t) on the level surface, the time rate of change of log(|V$|) on the 
level surface, and the surface divergence of flow velocity. 

The dynamic equations for global variables can be established once the transport 
equations for the scalar $(x,<) defining the level surfaces and the scalar f(x,t) are 
known. In the case of flames, the transport equation for a reactive scalar field may 
be written: 


Dt 


= Q<j> , 


where D/Dt = d/dt + u.V is the substantial derivative and 

1 d 


11 ^ 


f> dx 




W*(*) , 


(9) 


( 10 ) 


denotes the imbalance of diffusion and source terms. The transport equation for / 
is assumed to be in the generic form 


The iso-surfaces are defined implicitly as solution of 

= $(x,t) -y = 0 , 

and the unit normal vector of the iso-surface is defined by 

1 d$ 

~ |V$| dx Q ' 

The surfaces move in their normal direction, with the velocity 


( 11 ) 


( 12 ) 


(13) 


v a — v a -|- . 


(14) 
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The relative progression velocity V of the level surfaces is determined by the 
imbalance Eq. 10 of diffusion and source term in Eq. 9: 

1 d d$ 1 1 

V = • (15) 

p|V$|&r a v dx Q ; | VO j | VO | 

The dynamic equation for the global variable T* then follows in terms of volume 
integrals as 

tt - / w - * )|V41 [;£ (' A s£) + "■"'£;)] 

(16) 

This is one of several equivalent versions of the dynamic equation for the global 
variable in variable density turbulence. The equation for the surface area 

t) is obtained by choosing / = 1 (see Eq. 3) 

tt = / <**-*> H (fe - + (' r £) + ' *•) t: 

( 17) 

The corresponding equation in terms of surface integrals can easily be obtained 
from Eq. 17 by using Eq. 2. In particular, the surface stretch 


may be cast in the form: 


where 




£(!.+&) , 

Ii = J dA(x)4>,(x,t) , 

s.M 

I 2 = j dA(x)V(x,t)^ . 


In Eq. 20 4>, is a notation for the in-plane strain rate: 

, dv a dva 


( 21 ) 
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and the imbalance between reactive and diffusive effects f l <p has been written in 
terms of the relative progression velocity V by using Eq. 10. When the scalar $(.r , t) 
denotes the progress variable used in premixed combustion ($ = 0 in fresh gases 
and <£ = 1 in burnt gases), then by choosing as iso-surface the surface corresponding 
to the value of the progress variable used to locate the flame front, the well-known 
result is found that flame stretch is directly related to both contributions of strain 
rate acting along the flame front ( X \ ) and curvature linked with propagation of the 
premixed flame (I 2 ) (Candel and Poinsot 1990). 

Inspection of Eq. 17 shows that the surface area changes due to the strain rate in 
the tangential plane of the surface caused by the motion of the fluid and the dynamic 
change of the scalar field defining the level surface. Both phenomena may increase 
or decrease the surface area. The dynamics of the scalar variable defining the level 
surface appear in Eq. 17 as molecular diffusion and the source term. It is clear that 
both may be positive or negative, hence may increase or decrease the surface area 
as time evolves. The molecular diffusion term can be analyzed in more detail for 
incompressible flows with constant diffusivity T, since it is then proportional to 

1 _ dn a drip dn 1 n Q n^ d 2 $ 

|V$| dx a dx Q dx$ + dx 1 |V$| dx a dxp 

The first term on the right side is purely geometrical since it only depends on the 
mean curvature H 


1 dn a 

2 dx a ’ 

whereas the second term depends on the variation of the defining scalar normal to 
the level surface 


and it follows that 


5$ 

dn 


— n * V3> 


5 


= AH 2 - 2 

|V$| dx a |V4>| dn 2 

holds, which proves that the effect of molecular diffusion on the surface area is 
not definite: The first term on the right side always increases the surface area, if 
the mean curvature is nonzero, but the second term may decrease or increase it 
depending on the signs of the normal derivative and the mean curvature. 

2.3 Solution of the area equation for material surfaces 
For material surfaces Eq. 17 for the area of level surfaces is reduced to 


dA* 


(<M) 



dVg 

dXr 


— n Q np 


dvp\ 
dx a ) 


■S*(v>>0 


dt 


( 22 ) 
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with the initial condition 


A*(V,0)= J dA , (23) 

as the area of the surface $(£,0) = <p. Material surfaces enjoy the property that 
they are the one-to-one image of the initial surfaces. It is clear that Lagrangian 
variables are the proper choice for the formulation of their evolution. Denoting by 
r and a the Lagrangian variables defined as time and initial position of material 
points, the position 2L{ r ,a) at a later time is then the solution of 


= u a (r,X), ^o(0,a) = a Q , (24) 

where v a (t,x) is the Eulerian velocity and x are the usual Eulerian variables. The 
inverse transformation (from Eulerian field to Lagrangian field) a a = X” 1 (t,x) is 
unique and smooth as long as the Eulerian velocity field is smooth. The solution to 
Eq. 23-24 can be constructed if the co-area formula is applied 


A*{p,t) = J dx\V$\6($(x,t)-<p) , 

v 

and the fact that the level surface is embedded in a part of the flow field is used. 
Taking for T>(r) a subset of the flow field that contains the level surface at time r 
and using the fact that the mapping X_{r, a) is defined not only of the level surface 
but on the set D(r), we transform the volume integral to the Lagrangian frame as 
integral over X>(0) 

A*(v,t)= J </aJ(a, r)|V$|(a, t)£($Ql, t) — <p) , (25) 

D(0) 

where J denotes the Jacobian of the transformation given by 


„ . 1 dX a dXp dX y 

J(a,r) - 6 e a ^e St)u — — ^ 


(26) 


This Jacobian is the determinant of the Lagrangian deformation tensor J a p 


T _ 9X ° 

Jq/3 ~ dap ’ 

which is nonsingular for smooth flows with J > 0. The points on the level surface 
are materially invariant so 


£($(£, 0 -<?) = £($(«, o)-v0 . 


and it follows that the area at time r is given by 
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A*W,*)= j <ia|V*|(o,0)<(i(a,0)-v.)J(o,r)|||&^ 

V(Q) 

where da is the volume differential at time zero. Application of the co-area formula 
leads to a surface integral over the initial surface 

Mv,r)- J ^(0)-'(a,o}||g^ , (27) 

S*(v> °) 

where dA{ 0) denotes the surface differential at time zero. This is the desired solution 
of Eq. 23 and Eq. 24, as can be verified easily by differentiation. The solution Eq. 27 
can be recast in terms of the Lagrangian deformation tensor if the kinematic relation 
(Ottino 1989, section 2.7) 

dA(r) = (28) 

is applied. It follows at once that the solution Eq. 28 is thus given by 

A*(v,t)= f dA(0)J(a,T)fn a n^J~^J~f\ , (29) 

s*(v,o) 

valid for material surfaces. The version Eq. 29 of the solution shows the special 
property of material surfaces possessing nonsingular Jacobians. This is not the 
case for surfaces moving relative to the fluid, and Eq. 27 and Eq. 29 axe then not 
necessarily well defined. 

It is worth noting that for incompressible flows J = 1 holds and, assuming con- 
stant molecular diffusivity T, the result 

(30) 

s*(v>,0) 

follows from Eq. 27. It shows that the evolution of the area of material surfaces 
for incompressible flows is determined by the initial surface and the values of the 
scalar dissipation = TV$ • V$ on the surface. The integral can be computed if 
the Lagrangian scalar dissipation r) is known, which is often the hard part of 
the problem. The expression for the area Eq. 27 does not contain any statistical 
information as only surface integrals are involved. It is equal to the expected area for 
homogeneous turbulence but is different from expected values for non- homogeneous 
flows. 

2.4 Computation of the area of level surfaces 

Consider a scalar field $(x,*) > 0, such as enstrophy, temperature, density etc., 
that is a realization of a turbulent field. The equation of the iso-surfaces of $ is 
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z X 


FIGURE 1. Level surface of kinetic energy for a periodic viscous flow generated by 
the interaction of two vortex tubes, which are initially off-set and orthogonal. The 
flow is used as example for the computation of the surface area. 

then given implicitly by 54>(^) = {x £ V : $(x) = <£>}• It follows from Eq. 2 that 
surface integrals can be represented as volume integrals over a fixed domain T> (note 
that S$(<p) is not the boundary of the domain T>\ hence, the Gauss theorem is not 
applicable). The surface area is 

A*(^,f) = J dx\V$\8{$(x,t)-<p) (31) 

V 

in which the Dirac-pseudofunction selects the points in space where the level sur- 
face is located. This property indicates that an integral transform with respect to 
the range of the scalar <£ would remove the delta function. Since the scalar <l> is 
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q 

FIGURE 2. Surface area of level surfaces for the flow illustrated in Fig. 1. The sym- 
bols correspond to the triangulation method and the line to the Laplace transform 
method Eq. 33 and Eq. 34. 


nonnegative, it follows that the Laplace transform is the appropriate tool. Hence, 
the transformed surface area is then defined by 

oo 

^4*(s) = J difiAi(<fi)exp(-sif) , (32) 

0 

where s = u + iw is a complex variable. Using Eq. 31 we obtain 


which leads to 



A*(s) = j dx|V$|exp(— s$(x)) . (33) 

v 

This is a powerful result with both practical and theoretical implications: the 
Laplace transform of the surface area of level surfaces for all level values can be 
computed without computing the location of a single level surface. Furthermore, 
this computation is much more efficient than methods based on triangulation of 
level surfaces that keep track of the connectedness of the triangles. It is also useful 
for theoretical considerations since the integrand in Eq. 33 is a well defined function 
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for finite s and not a distribution, and is accessible to the methods of classical 
analysis. Once the Laplace transform has been determined, the inverse transform 
produces the desired surface area according to 

w + ioo 

• 4 *( v >)= 2 ^ J dsA$(s) exp(stp) . (34) 

u — ioo 

The numerical computation of the integrals in Eq. 33 and Eq. 34 requires some 
care. It is easy to see that the range of w in the independent variable s = u + iw in 
Eq. 33 and Eq. 34 is limited by the highest wave-number that can be resolved on 
the lattice of discrete values of <p 


u> max ||A$(r,t)|| < 7T , 

where || * || is an appropriate norm and A$ denotes the difference of $ between two 
neighboring lattice points. Choosing w max significantly larger than this limit leads 
to integrals over rapidly oscillating functions which diverge. For the present case 
the maximum norm was initially chosen, which turned out to be too restrictive, and 
the limit was relaxed to twice the value given above. This aspect of the transform 
method requires further analysis. 

As an example, the scalar kinetic energy $(x,f) = q(x, t) = v Q v Q generated as a 
numerical solution of the Navier-Stokes equations for the interaction of two initially 
orthogonally offset vortex tubes (Re = 1500) at time t = 1.375 is used to illustrate 
the transform method outlined above. One of the level surfaces is shown in Fig. 1 for 
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the level value which is 48.78% of the maximum. It has the Euler number (Kollmann 
and Chen 1994) \ — 10, which corresponds to the genus </ = 1 - jx = 6 °r a 

sphere with six handles, and the fragmentation N = 1, which means that there is 
only one surface component. The intersections of the level surface (the flow field 
is periodic in all three directions) with the cubic flow domain are artificially closed 
with lids. The area of the lids is not included in the results for the surface area. 
Two methods for the computation of the surface area are applied: Triangulation of 
the level surfaces, which determines the connectedness information for all triangles 
and thus allows also the computation of the topological properties x and N, and the 
Laplace transform method given by Eq. 33 and Eq. 34. The results in Fig. 2 (linear 
scale, the full line is the transform method and the symbols the triangulation) and 
Fig. 3 (semi-logarithmic plot) show that good agreement between the two methods 
is achieved. The advantage of the Laplace transform method becomes apparent if 
the CPU times are compared: The Laplace transform method (using 100 points 
on the ^-axis) is about 1000 times faster than the triangulation method (using 40 
points). 

2.5 Surface density function 

In modeling purposes one wishes to seek out closed equations for point quantities 
that can be easily discretized and introduced in fluid mechanics codes. The density 
of iso- surface is introduced in this section as a quantity that characterizes, in mean, 
the behavior of iso-concentration surfaces for a reactive scalar field $(r ) that follows 
Eq. 9. 

The basic relation (Eq. 2) set out in Section 2.1 suggests that for / = 1 a local 
and instantaneous measure of the density per unit volume of iso-surface is provided 
by 


S(^;^0=|V$U)|^(x)-^) . 


Indeed from Eq. 2 



(35) 


and £(<r>;x, 0 corresponds to the ratio in the limit 6V — ► 0. This expression 
indicates that, within a turbulent environment, the mean surface density function 
(sdf) E(<^;x, t) is the product of the expected value for the magnitude of the gra- 
dient of $(x, t), conditioned on being on one iso-surface $(x, f) = ip, times the 
probability of being on that surfa ce. Indeed introducing the pdf of $(x, t) = 
namely P{p>\x,t) = 6(p - $(x,t)), one may write: 


S(^;x,<) = |V$(x,f)| %-*(*,<)) = 


( 


| V*(x,t) I $(x,0 
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where V$(z, /) | |$(x, t) = pj denotes a conditional average. (This definition of 

surface density has been previously used by Trouve and Poinsot (1994) to evaluate 
flame surface density from DNS.) 

All information contained in the sdf is also available through the joint pdf of $ 
and its gradient; models have been proposed for this joint pdf in the case of reacting 
and non-reacting flows (Meyers and O’Brien 1981, Dopazo 1994, Fox 1994). The 
iso-surface is one of the standard concepts used to investigate complex phenomena 
in premixed turbulent flames (Veynante et al. 1994); therefore, despite the detailed 
statistical description involved in this joint pdf, it is interesting to restrict the 
present study to the specific equation for E(y?;x,t). 

An exact transport equation for the surface density function S(^;x, t) may be 
written (Vervisch et al 1994): 


= ((-Hn:Vu) | V$Gr,f) | |<K*,0 = ^ P(^x,t) 

— (n.Vfi$(x,f)|$(:r, t) = P(<p',x,t) 

- ^ (ft* (x, t) $(x, t) = <pj P(p; x, t) 


(37) 


It is important to note that in Eq. 37 p is an independent variable (sample space 
variable), and Eq. 37 is valid for all the iso-concentration surfaces defined from the 
progress variable field. When the various terms contributing to the sdf evolution are 
transposed to a surface formalism, they are equivalent to those of Eq. 4 or Eq. 17. 

2.6 Unified treatment of pdf and sdf equations 

Pdf modeling has been shown to be an efficient tool for the computation of turbu- 
lent non-premixed flames (Kollmann 1990, Borghi et al. 1991). The pdf approach 
has a great potential because it allows for including the description of the chemistry 
in a closed form, the modeling issue related to the evaluation of the mean burning 
rate being reduced to the estimation of the small-scale mixing. However, in the 
case of turbulent premixed flames, classical mixing models may fail because of the 
particular properties of the reactive and diffusive zones. Indeed, mixing is a pre- 
requisite for reactions in non-premixed systems as is preheating in premixed flows; 
therefore, the resulting predictions depend crucially on the properties of the mixing 
model, no matter how accurately the chemical reactions can be represented. 

The pdf equation may be written (Kollmann 1990): 


£>P(y;x,Q 

Dt 


d_ 

dip 


(ft+(x,*) $(x,<) = yj) P(<p;x,t) , (38) 


where the bar over the substantial derivative indicates that turbulent transport 
effects are included. By introducing the dissipation of the scalar field = T | | 2 , 
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the modification of the pdf in composition space may be written in an alternative 
form 


£>P( ^;- ,t) = ~ f = <p) P(v;x,t) 

Dt *£-1 > , J (39) 

- [ (' e * | ' » 0 = *p ) p (w * ) 

In Eq. 39 the chemical source is closed, but the mixing term involving second order 
derivative in composition space is unknown when a one-point description is adopted 
for the turbulent flow. This term is actually representative of the distribution of 
the magnitude of gradient of the scalar field along the iso-concentration surfaces. 
It is therefore natural to look for a closure of this term that takes advantage of 
the information included in E(^>; x, t) concerning the dynamics of these iso-surfaces. 
This approach leads to the development of a PDF-SDF model in which a system of 
equations consisting of the transport equations for the pdf and the sdf is solved. 

The sdf equation can be recast in a form that emphasizes the similarity to the pdf 
equation if a relation between conditional and surface means is established. Surface 
means may be defined by 





) 



J dA(x)Q$(x,t) , (40) 


where Eq. 2 has been used when E (<p\x,t) =| V$ | 6(<p — $) is the instantaneous 
surface density function. 

The conditional mean of may be written: 


$(*,<) = <?) = 


t)6($(x, t) — p)) 
6(${Xyt) ~ if) 


(41) 


A relation to the surface mean can be established for homogeneous flows if the 
co-area formula Eq. 2 is applied to Eq. 41. It follows for homogeneous flows that 
the conditional mean is the ratio of volume integrals which become surface integrals 
using Eq. 2 


= 9) = -^-7 

V J I dA(x) 


I V*| 

T 
I V*| 


which is recognized as the desired relation between conditional and surface means 
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(i}.j>(x,<) $(x,<) = = 


(n^if 2 •**=*) 

(iWT 


Both sides are defined for non-homogeneous flows but are not equal since the 
right side is the ratio of geometrical means Eq. 40 and not statistical expectations. 
Conversly, the surface mean can be expressed in terms of conditional means 


A<t>=vj == 


(ft,j.(x,f) I V$ I $(x,t) = y>) 

(l V* | | *(*,*) = ¥>) 


(42) 


valid for homogeneous flows. 

The application of the previous results to homogeneous flows converts the sdf 
equation (Eq. 37) in the form 


as (ip-,x,t) 

dt 


d_ 


(ft*|.4* =v ,) S(y>;x,f) + (^»|v4*= v ^ S(v?;x,<) , (43) 


where d a denotes the total rate of generation of surface area per unit volume 

8 a = —nn:Vv — . 

Comparison of the pdf Eq. 38 with the Eq. 43 for the mean sdf shows that the 
sdf is transported in scalar space by the same mechanism as the pdf, but is also 
created or destroyed by the tangential strain rate and by the flux of the imbalance 
between diffusion and reactions normal to the level surfaces. 

The equations for the pdf and the sdf in homogeneous flows pose two closure 
problems: the fluxes in scalar space and the surface mean of the total generation 
rate of the sdf. The particular structure of premixed flames and the previous analysis 
suggests dealing with mixing and chemical reactions in a unified way by developing 
a closure model for their combined effects. Pdf and sdf are moved in scalar space 
by the conditional mean and the surface mean of the same quantity, namely the 
imbalance between diffusion and reactions. 


3. The internal structure of premixed turbulent flames (DNS) 

A database of turbulent premixed flames has been generated by using a fully 
compressible direct numerical simulation code (Guichard, Vervisch and Domingo 
1994). These simulations are three dimensional, and fields are computed by using 
a 129x129x129 grid. The numerical approach chosen has been proposed by Lele 
(1990), Poinsot and Lele (1992). 
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FIGURE 4. Snapshots of the flame surface contour = 0.8). 


In previous studies of premixed and non-premixed turbulent flames at CTR 
(Trouve and Poinsot 1994, Vervisch 1992), initial conditions were generated by 
using a laminar and planar flame on which a homogeneous turbulent field was su- 
perimposed. As a consequence, prior to being distorted by the turbulence, the 
structure of the reactive zone corresponded to the one of a laminar flame. Globally, 
the computational domain was divided into fresh and burnt gases separated by a 
propagating turbulent flame. These simulations have been proven to be an efficient 
tool for studying, through the flame brush, the topology of one iso-surface identified 
as the flame-surface and corresponding to the location of intense reactive activity. 

In order to increase the number of possibilities of sampling for all the values of the 
progress variable field, in the present simulations the initial scalar field is prescribed 
according to a pdf and to a turbulent spectrum. It is then possible to generate 
as initial condition a homogeneous progress variable field with pockets of burnt 
gases and fresh gases. A turbulence velocity field is added to this scalar field, and 
turbulence is allowed to decay. These pockets axe convected by the turbulence, and 
reactive layers are created. To avoid spurious effects that can emerge from this initial 
condition, the procedure used to initialize the scalar field insures that the transition 
between fresh and burnt gases is achieved according to a smooth distribution. Non- 
reflecting boundary conditions are used on each side of the computational domain, 
and pressure waves can leave this domain. 
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Surface means in flame units 



FIGURE 5. Distribution in progress variable space of surface means, 
turbulent flame; ; laminar flame. 


From the initial condition, the progress variable field is convected by the turbu- 
lence, and propagation of premixed flames is observed. A single-step chemistry with 
a reaction rate cast in the form of an Arrhenius law represents the chemical activ- 
ity. The ratio of temperature between burnt and fresh gases is five, the Zeldovich 
number has a value of eight. The initial ratio between the rms turbulent veloc- 
ity and laminar flame speed is the order of three, and the value of the Damkohler 
number based on the Kolmogorov time scale and the laminar flame time is about 
unity. The initial turbulent Reynolds number based on the integral length scale is 
about seventeen, while the Reynolds number based on the initial Taylor microscale 
is about fifteen. The formalism which has been developed above is used to analyze 
data from this calculation. 

After one flame time, the flame surface exhibits a very complex and convoluted 
shape that can be observed in Fig. 4. 

Fig. 5 presents the distribution in progress variable space of surface means as 
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defined in Eq. 40 for various quantities; their corresponding distribution across a 
one- dimensional steady unstrained laminar flame is also plotted. The distribution 
of surface means for the diffusive budget of the progress variable ^(V.(prV$(x, t))), 

of the chemical source and of the magnitude of the gradient | V$(x, t) | follow 
the laminar response except in the preheat zone ($ < 0.5), where, because of the 
turbulence activity, the flame is not able to sustain the flamelet structure. The 
turbulence clearly affects the progress variable gradient in the preheating region. 
As a consequence, the surface mean of the iso-surface relative progression velocity 

drops in the preheat zone and in the burnt gases. One may observe with the help 
of these plots how the imbalance between chemical source and diffusion leads to 
the propagation of the various iso- surfaces. However, even though the surface mean 
value of this velocity is close to the laminar flamelet value, the scatter plot (Fig. 6) 
of this quantity shows that locally the effect of the turbulence can be very large in 
this DNS. This result suggests that there is not a unique deterministic relationship 
between progress variable and propagation speed. The same observation has been 
made in connection with the imbalance between reaction and diffusion as well as 
for the magnitude of the gradient. 

Fig. 6 also displays the pdf of the progress variable at time t = 0 and the pdf 
after one flame time; a bimodal shape is observed in accordance with the existence 
of the flamelets seen in Fig. 5. The effects of the turbulence on the preheating zone 
is emphasized in Fig. 6 where the area of the iso-surface is increased for small value 
of the progress variable. 

All of these observations suggest the existence of a flamelet zone within the tur- 
bulent flame brush, where mean conditional surface values follow the laminar flame 
structure. However, the existence of a preheating region is also observed, and the 
properties of the mixing in this zone may control the behavior of the entire flame 
structure. 

The total flame stretch, defined for the iso-surface corresponding to the progress 
variable that gives the peak value of the reactive activity, is at a time in the sim- 
ulation corresponding to Fig. 5 and Fig. 6 positive when the mean value for the 
progress variable evaluated over all the computational domain is $ = 0.55. Flame 
stretch is known from DNS or experimental observation (Trouve and Poinsot (1994) 
Veynante et ai 1994) to evolve from positive to negative values as a function of the 
mean progress variable. In accordance with these observations, later in time, flame 
stretch becomes negative when the mean value of the progress variable within the 
domain becomes larger. 

4. PDF-SDF modeling 

DNS observation suggests that the premixed turbulent flame may be organized 
in progress variable space in three zones: a preheat zone and a burnt gas zone 
where mixing needs to be treated separately from the slow chemical reaction, and 
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Flame unite 



FIGURE 6. Clockwise from upper left: Pdf of the progress variable ( , pdf; 

, initial pdf), scatter plot of iso-surface relative progression velocity, area of 

iso- surfaces, and iso- surface density function. 


a flamelet zone where the mean conditional surface value of the velocity of the iso- 
surface meets the one of a laminar flame and where laminar flamelet behavior may 
be invoked. Because of the necessity of introducing the chemical rate within both 
the preheat and burnt gases zones, Eq. 39 is appropriated for the pdf, while Eq. 38 
is well suited for the flamelet zone where the convection in composition space will 
be deduced from flamelet behavior and represents the laminar flamelet 

isosurface density. It is zero when no laminar flamelet burning zone exists. 

A direct approach to close the pdf equation can be developed by prescribing for 
the imbalance the value corresponding to the planar unstrained laminar flame 
. When combining with a classical mixing model, a closed pdf mixing term 
emerges (Anand and Pope 1987). Here, our objective is to take advantage of both 
pdf and sdf transport equations at the same time. The value of corresponding 
to the laminar flamelet structure will be used, but in conjunction with strain-rate 
curvature and related mixing effects in both pdf and sdf equations. 
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Flame units 



Progress variable 

FIGURE 7. Comparison between modeled sdf terms (valid within the flamelet 
range) and DNS. : DNS; : proposed model. 

4-1 Modeling for the sdf equation 

In this work, the turbulent convection will not be considered and its investigation 
may be found in this book elsewhere (Trouve et al 1994). 

From the above analysis, the equation for the flamelet sdf 

%>;*,<) = (| V* I |*(x,<) = ¥>) P(<p;x,i) , (44) 

may be written: 


Term 1 
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FIGURE 8. Comparison between modeled pdf terms (valid within the flamelet 

range), top, and DNS (flame units), bottom. Top: , conditional mean value; 

, conditional surface mean value. Bottom: , DNS; , proposed 

model. 


— *)|$(x, t) = (p'j P(y?;x,f) 


(45) 


Term 2 


_a_ 

dip 


(ft$(x,i) I V$(x,t) I $(x,i) = <r>^ P(v,x,t) 


Term 3 


The strain rate term 


(< <t>» is 


is representative of a physical process that 

affects in a similar way all the iso-surfaces of the scalar field. As presented in 
Fig. 7 (Term 1) this contribution does not strongly depend on the progress variable. 
Therefore, it is possible to speculate that model closures developed for the in-plane 
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strain rate in the framework of the flame-surface density approach can be extended 
to the sdf equation (Meneveau and Poinsot 1991). Attention will be especially 
focused on Term 2 and Term 3. 

For Term 2 one may prove that the expression 


- <)|$(x, t) = <pj , 

does not explicitly contain curvature effects (Vervisch et al. 1994). The contribution 
of curvature of the field is actually present through the pdf P(<^>; x, t) that multiplies 
this term in the sdf budget. Hence it is possible to build a closure by using the 
flamelet structure. If S/f(<^>) denotes the known laminar flame speed, evaluated in 
the premixed laminar unstrained flame as 


Sfiif) = 


ft*, 

I Vy> I" G( V ) 


where G((p) is the known magnitude of the gradient of the progress variable within 
the laminar flame, then one may write: 


- ^n.VfU(r,t)|$(r,0 = V 7 ) P(wz,t) - ^(v 7 ! 2L,*)G(<p)^(Sfi(<p)G(<p)) . 

Evaluated from DNS the exact expression (left side) and the proposed model 
(right side) are displayed in Fig. 7 (Term 2). Good agreement is found according 
to a possible definition of the flamelet range [po,<pi] = [0.4, 0.8]. 

The same approach is followed for Term 3, and this leads to 




PiWZyt) 




The validity of this modeled expression within the flamelet range is evaluated 
in Fig. 7 (Term 3). One may note from Fig. 5 that according to the distribution 
of the progress variable gradient, flamelet range could be [v^o^i] = [0.6,1] while 
the iso-surface relative progression velocity follows flamelet behavior for the range 
[0.4, 0.8]. The same trend is observed in Fig. 7 where the modeled term related to 
the propagation (Term 2) follows its exact DNS expression in the range [0.4, 0.8], 
while good agreement is obtained for the closure related to mixing in the range 
[0.6, 0.8]. 

The closed transport equation for the sdf is only valid within the flamelet zone; 
however, information coming from the preheat and burnt gases zones is present 
through the pdf P(y?;£, <). The partially closed equation contains one source term 
and one term representing transport in composition space and reads: 
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+ v.(Sfo£,«)v) = 

( 1 O ^ \ 

-P(^;x,*)5 /< (^)^G !2 (v>) + 5/i(v>)G 2 (¥>)^?(¥);i,i)J > 

(46) 

where the turbulent transport effects are included in the second term on the left 
hand-side. The function W(v?o? V^i) insures that when flamelets are not present the 
RHS of Eq. 46 goes to zero (W(v?o J V ? i) equal to zero when the flamelet range 
[v^Oi^i] is not defined and unity when flamelets are present). 

4-2 Modeling for the pdf equation 

To treat both flamelet part and mixing zones, from Eq. 38 and Eq. 39 the exact 
transport equation for the pdf may be cast in the form 


—^f 1 — - -<*(¥>; $(*><) = 9^ P(v;x>*) 

- ( 1 . - 7(ip;x,t) 

~ (1- - j ^«*|*(*><) = £><) » 


(47) 


where a(^;£,t) is a smooth function equal to unity in the flamelet zone of p space 
and to zero in the outer flamelet regions. A possible shape for the function a is 
proposed in Fig. 7. The way in which the flamelet range ls *° be 

determined will be considered later. 

Within the out-of-flamelet zone, a classical mixing modeling (Kollmann 1990) 
can be used with complex chemistry to evaluate the contribution of the reactive 
activity. 

In the flamelet zone, from Fig. 8 it is observed that the surface mean of the 
imbalance ft* is almost equal to its conditional mean; from this we deduce the 
following modeled expression for the flamelet part of the pdf equation 


H *<*•*)-*) - ( n ‘M - - 


and hence 


~k [H* - *) “ - (ffirij) 


9 a (WiiO 

■x- 5/|(v>)G (*>) -=- — 

ay S(v;*.0 


( 48 ) 
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is introduced as a closure in Eq. 47. The validity of this expression is evaluated 
against the DNS results in Fig. 8. 

The determination of this flamelet range is still unproven at this stage. However, 
we speculate that a criterion based on the imbalance between the reactive- mixing 
terms and the flamelet term (Eq. 48), both present in the pdf equation, may emerge 
to determine dynamically the values and 

The partially closed system may be written: 


- My* o,<^i) + S f i(v)G 2 (v)^P(wx,t)j 


•P-PQ;g,0 

Dt 






Sj\W)G (<*>)=; 


E(v?;x,<) 

- (1- - a(v?;x,f))^- $(*,<) = ^ P(<p\x,t) 

- (l. - a(<p-,x,t))-^- . 


in which closed expressions for the strain rate and for turbulent transport remain to 
be included. The LMSE (or IEM) mixing model (Dopazo 1994) has been chosen for 
the outer flamelet zone and r t is a characteristic mixing time scale. Monte, Carlo 
simulation is appropriated for the pdf equation (Pope 1985), while a continuous 
solution can be considered for a determinate numbers of bins for which the sdf 
equation is solved. 


5* Conclusion 

Because the motion of iso-surfaces is strongly connected with the transport phe- 
nomena at all scales of the turbulent flow, the behavior of scalar iso-surfaces in 
multi-component reacting mixture is a topic of fundamental interest. It is possible 
to develop some mathematical tools that describe the properties of these surfaces. 
In particular, a transport equation for global variable has been established and a 
solution for material surfaces evolution has been proposed. When introducing a 
surface density function, a possible link between the well-known pdf approach for 
turbulent combustion modeling and surface modeling emerges. From this point an 
attempt is made to develop a unified treatment of pdf and sdf for premixed turbu- 
lent flames. Direct numerical simulations results are used to analyze the internal 
structure of turbulent premixed flames and to estimate the validity of the proposed 
modeling. 


Ack nowledge me nt s 

The authors would like to thank their CTR hosts G. Ruetsch and J.-M. Samaniego 
for helpful assistance. J. Reveillon and L. Guichard are thanked for their help in the 



PDF-SDF modeling for premixed turbulent combustion 


149 


development of the DNS database. Authors have benefited from many discussions 

with the participants of the “combustion group” during the CTR summer program. 

REFERENCES 

An AND, M. S., POPE, S. B. 1987 Calculations of premixed turbulent flames by 
pdf methods. Comb, and Flame. 67, 127-142. 

BORGHI, R. 1988 Turbulent combustion modeling. Progr. Energy Comb. Sci.. 14 , 

BORGHI, R., VERVISCH, L., Garr£ton D. 1991 The calculation of local fluctua- 
tions in non-premixed turbulent flames. (Eds. Carvalho M. G ., Lockwood , F., 
Taine , J.), Springer Verlag , 88-118. 

BRAY, K. N. C. 1990 Studies of the turbulent burning velocity. Proc. R. Soc, 
Lond. A, 431, 315-335. 

CANDEL, S. M., POINSOT, T. J., 1990 Flame stretch and the balance equation 
for the flame area. Combust Sci. and Tech . 70, 1-15. 

CONSTANTIN, P. 1992 Remarks on the Navier-Stokes equations. To appear. 

DOPAZO, C. 1994 Recent developments in pdf methods. In Turbulent Reacting 
Flows (Eds P. A. Libby and F. A. Williams), Academic Press London, 375-474. 

Fox, R. O. 1994 Improved Fokker-Planck model for the joint scalar, scalar gradient 
pdf. Phys. of Fluid. 6 (1), 334-348. 

GuiCHARD, L., Vervisch, L., Domingo, P. 1994 Numerical study of the inter- 
action between a mixing zone and a pressure discontinuity. AIA A 95-0877 , to 
be published. 

KOLLMANN, W. 1990 The pdf approach to turbulent flow. Theoret. Comput. Fluid 
Dynamics. 1 , 249-285. 

KOLLMANN, W., CHEN, J. H. 1994 Dynamics of the flame surface area in turbulent 
non-premixed combustion. Twenty-Fifth International Symp. on Comb. 

Lele, S. K. 1990 Compact finite difference schemes with spectral like resolution. 
Center for Turbulence Research, Stanford Univ./NASA Ames, Manuscript 107. 

Meneveau, C., POINSOT, T, 1991 Stretching and quenching of flamelets in pre- 
mixed turbulent combustion. Comb, and Flame. 86, 3311-332. 

Meyers R. E., O’Brien, E. E. 1981 The joint Pdf of scalar and its gradient at 
a point in a turbulent fluid. Comb. Sci and Tech . 26, 123-134. 

Marble, F. E., BROADWELL, J. E. 1977 The coherent flame model for turbulent 
chemical reactions. Project Squid Report, TRW-9-PU, TRW, El Secundo. 

Maz’ja, V. G. 1985 Sobolev spaces. Springer V. Berlin . 

OTTINO, J. M. 1989 The kinematics of mixing: Stretching, chaos and transport. 
Cambridge Texts in Applied Mathematics, Cambridge University Press 1989. 

PoCHEAU, A. 1992 Front propagation in a turbulent medium. Europhys. Lett. 20 
(5)„ 401-406. 



150 


L. Vervisch, W. Kollmann, K . N. C . Bray , & T. Mantel 


PoiNSOT, T., Lele, S. K. 1992 Boundary conditions for direct simulations of 
compressible viscous flows. /. Comput. Phys. 101(1), July 92. 

Pope, S. B. 1985 Pdf method for turbulent reacting flows. Comb . Sci. and Tech. 
11, 119-195. 

TROUV6, A., PoiNSOT, T. 1994 The evolution for the flame surface density in 
turbulent premixed combustion. J. Fluid Mech. To be published. 

Trouv£, A., Veynante, D m Bray, K. N. C., Mantel, T. 1994 The coupling 
between flame surface dynamics and species mass conservation in premixed 
turbulent combustion. In Proceedings of the 1994 CTR Summer Program , 
Center for Turbulence Research, NASA Ames/Stanford Univ, 

VEYNANTE, D., Duclos, J. M m PlANA, J. 1994 Experimental analysis of flamelet 
models for premixed combustion. Twenty-Fifth International Symp. on Comb . 

Vervisch, L, Bidaux, E., Bray, K. N. C., &, Kollmann, W. 1994 The 
transport equation for an iso-surface density function. Submitted to Phys . of 
Fluids . 

VERVISCH, L. 1993 Study and modeling of finite-rate chemistry effects in turbulent 
non-premixed flames. Annual Research Briefs 1992 . CTR, Stanford U./NASA 
Ames. 



